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Abstract. - We investigate the steady state heat current in two and three dimensional isotopically 
disordered harmonic lattices. Using localization theory as well as kinetic theory we estimate the 
system size dependence of the current. These estimates are compared with numerical results 
obtained using an exact formula for the current given in terms of a phonon transmission function, 
as well as by direct nonequilibrium simulations. We find that heat conduction by high-frequency 
modes is suppressed by localization while low-frequency modes are strongly affected by boundary 
conditions. Our heuristic arguments show that Fourier's law is valid in a three dimensional 
disordered solid except for special boundary conditions. We also study the pinned case relevant to 
localization in quantum systems and often used as a model system to study the validity of Fourier's 
law. Here we provide the first numerical verification of Fourier's law in three dimensions. In the 
two dimensional pinned case we find that localization of phonon modes leads to a heat insulator. 
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Transport in random media is a topic of great current 
interest flH^. Here we study the effect of isotopic mass 
disorder and boundary conditions on heat conduction and 
localization in two (2D) and three (3D) dimensional har- 
monic crystals, where the properties of the ordered system 
are known exactly. These are systems of N d (d = 2, 3) 
atoms in contact, at their surfaces perpendicular to the 
x-axis, with heat reservoirs at different temperatures. We 
shall focus mainly on the dependence of the heat flux J 
on N. When Fourier's law holds then J ~ TV -1 but this 
is known to be violated in computer simulations [8] and 
some experiments (9 10 in ID and 2D systems where 
one finds J ~ N~^ with /i / 1. Earlier work [7,11] on 
heat conduction in glassy harmonic systems computed the 
frequency-dependent thermal diffusivity using the Green- 
Kubo formula for different system sizes. However their 
study did not directly address the question of asymptotic 
size dependence of conductivity. 

We will first describe the precise model studied, then 
give our heuristic arguments for calculating the N depen- 
dence of J and finally present the numerical results. 

We consider simple cubic lattices with displacements at 



each lattice site n (n„ = 1, 2, ...,N for v = 1,2, ..., d) given 
by a scalar variable x n . In the harmonic approximation 
the system Hamiltonian is given by 



iV-l 



" - E?-n+EE^-W 



ni — 1 n' e 



E7fc+4,,)+E^n (1) 



where e refers to the d unit vectors and we have denoted 
n = (?ii,n'). We impose periodic boundary conditions 
(BCs) in the directions v = 2, 3, ...,d and consider two 
different BCs in the direction of heat conduction (y = 1): 
(i) fixed BCs k' > and (ii) free BCs ft' = 0. Pinning 
refers to the case with k > 0. We consider binary mass 
disorder with equal number of particles of masses fh — A 
and fh + A distributed randomly on the lattice sites. The 
equations of motion of particles in the bulk ( 1 < n\ < N 



) are given by: 



2^, ±iji k(x n 



k, n 



The particles at the surfaces m = 1 and n\ — N are 
connected to heat reservoirs, at temperatures Tl and Tr 
respectively. These are modeled by white noise Langevin 
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equations. This means that particles at ri\ — 1 and n% = 
N, have additional forces given respectively by (— k'x n — 



7-in + Vn) an d (—k' 



7^n + ^n)i where rj and rj 



are noise terms with strengths proportional to Tj, and Tr 
respectively and to the friction constant 7. If we model 
the heat reservoirs themselves by infinite ordered harmonic 
crystals then Langevin type equations for the system 12 
are obtained on eliminating the bath degrees of freedom. 
The two different BCs then emerge naturally. Fixed BCs 
correspond to reservoirs with properties different from the 
system (e.g. different spring constants) while free BCs 
correspond to the case where the reservoir is simply an 
extension of the system (without disorder) 6, 13 14 . We 
note that fixed BCs are more realistic. 

In the noncquilibrium steady state, the heat current per 
unit area (assuming unit lattice spacing length) from the 
left to the right reservoir is given by 



12 



15 



J 



AT 
2nN d ~ 1 



dujT(uj) , 



(2) 



where AT = Tj, — Tr and T(to) is the transmission co- 
efficient of phonons at frequency w from the left to the 
right reservoir. Our interest here are in the disorder av- 
eraged transmission T(ui) — \T{uj)}/N d ~ l and the current 
J=[J]. 

The phonon and electron localization problems are 
closely related. For the system without baths consider the 
displacement field a n (p) for the p th normal mode. This 
satisfies the equation: m n uj 2 a n = (2d + k )a n — J2e a n+e- 

By introducing variables ipn(p) = m n a n(p)i v a = {2d + 
^o)/ TO n and t n ,i = l/(m n rni) 1 / 2 for nearest neighbour 
sites n, 1 this equation transforms to a Schrodinger-type 
equation uj 2 ip n (p) = v n i[) n (p) - X^Vi^Mp) ■ The near- 
est neighbour hopping t n ,\ and on-site energies v n are now 
correlated random variables. For the unpinned case with 
k = 0, translation invariance gives rise to extended modes 
at low frequencies. The pinned case has no translational 
invariance and is closer to the usual electron localization 
problem. 

The effect of disorder on heat conduction in a harmonic 
crystal manifests itself in two ways: (i) Anderson local- 
ization [16] of phonon modes of frequency w will make 
them non-conducting, (ii) non-localized phonons will be 
scattered by the impurities. From Eq. pi) the net current 
is given by the integrated transmission and for large sys- 
tem sizes the integral is over the range of normal modes. 
It is natural to classify normal modes as localized, diffu- 
sive or ballistic. For modes localized on a length scale £, 
T(oj) ~ e~ N ' £ . This £ depends on the phonon frequency 
and low frequency modes for which £(uj) ~ N will therefore 
be carriers of the heat current. Diffusive modes are spa- 
tially extended but non-periodic and T(uS) ~ 1/N. Bal- 
listic modes are extended and approximately periodic and 
their transmission is TV- independent. 

A renormalization group study in a disordered contin- 



2D all non-zero frequency phonons are localized. They 
studied the spreading of an energy pulse to define a fre- 
quency dependent diffusivity Dq(oj), From the behaviour 
of Dq(uj) under renormalization one can obtain a differ- 
ential recursion relation for the resistivity. This relation 
shows that in the large system size limit the RG flow is 
towards infinity for dimensions < 2 in contrast to d = 3 
where the flow is towards zero as long as u> is less than 
some fixed value, independent of system size. Hence all 
finite-frequency modes in one and two dimensions are lo- 
calized. From the differential recursion relations one finds 
that for d = 1 and d = 2 the localization length in u — > 
limit diverges as ~ 1/oj 2 and ~ e 1 ^ respectively. In 3D 
there exists a frequency, independent of system size, above 
which all states are localized while states below that fre- 
quency are extended. Hence for a system of size N there 
will be a cut-off frequency uj\\ (which depends on N for 
d < 2) above which all the modes are localized. In differ- 
ent dimensions cu^ is given by 

lu^ ~ N~ 1/2 for d = 1 

- [log(N)]- 1/2 for d = 2 (3) 

~ nonzero value independent of N for d = 3 

For extended modes with w < w^ we use kinetic theory 
to determine the cut-off frequency ui^f below which states 
are ballistic. Rayleigh scattering of phonons gives a mean 
free path £k(u) ~ w~( rf+1 ) in d-dimensions. This gives 
w f = AT-i/(d+i) below which £ k (uj) ~ N. 

The net current in the system consists of both ballistic 
[Jb ~ Jq c du)T(uj)\ and diffusive [Jjj ~ f^ dujT(uj)) con- 
tributions, J = Jb + Jd- One other crucial observation, 
consistent with the numerics, is that T(uj) in the ballistic 
region has the same form as for an ordered system, which 
is sensitive to boundary conditions. 

In II? both uj^,uj^ ~ N^ 1 / 2 and numerical studies 
show that all extended states are ballistic [l8J. For free 

TV -1 / 2 

BC T(uj) — > const, as u — ¥ and hence J ^ J a dw ~ 

JV-V2. For fixed BC T(u) - uj 2 leading to J - A~ 3 / 2 . 
The dependence of the form of T{uS) on pinning was dis- 
cussed in 



19 in detail. Thus surprisingly one finds 



strong dependence on BCs. In the presence of pinning, the 



low frequency modes are removed and J ~ e 



-cN 



These 



results agree with earlier rigorous and numerical work on 
this system 13 15 18 20 . We now calculate the asymp- 



uum elastic model by John etal 17 found that in ID and 



totic system size dependence in 2D and 3D. 

Unpinned lattice: To find the ballistic contribution we 
note that for the ordered lattice T(lS) ~ uj d ^ 1 as u> — >• 
for free BCs and as T(u>) ~ cu d+1 for fixed BC. This can be 
seen in the following way. We write n = (m,n') and q = 
(<72,<?3, —,Qd) with q a = 2j ^ L where n goes from 1 to N. 
Now if we define x ni {q) = JV^- 1 )/ 2 £ n , x (ni , n -) e^-"', 
then one can show that, for each q, Xi 1 (q) satisfies a 
Langevin equation corresponding to a one-dimensional 
Hamiltonian with the onsite spring constant k replaced 
by -Mq) = fco + 2k(d — 1 — J2 a =2 d cos (la))- This means 
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Fig. 1: Plot of the A?-dependence of disorder averaged current 
J. (a) 2D: Inset shows result for pinned case. For free BC 
A = 0.8 and for fixed BC A = 0.95. For pinned case A = 
0.4, k = 10.0. (b) 3D: For all cases A = 0.8 and for pinned 
case k = 10.0. Error bars show standard deviations due to 
sample-to-sample fluctuations and are very small except in the 
2D pinned case. 



that the problem of heat conduction in a d-dimensional or- 
dered harmonic lattice can be related to heat conduction 
across N d ^ 1 independent ordered harmonic chains with 



different onsite potentials 21-23 . Hence the transmission 
coefficient T(u) for the d-dimcnsional lattice can be ex- 
pressed as a sum of the transmission coefficients of the ID 
ordered chains. Using this result and the analytic form of 
the transmission coefficient for the ID chain 
in the N — > oo and w — >• limit 



23 we find 



T(w) 



u! for open boundary condition 

U) + for fixed boundary condition. 



Hence we get the ballistic contribution to the total current 

K 

density (for the unpinned case) as: Jb ~ J " dio to d ~ l ~ 
N -d/(d+i) for free BC and Jb _ J£° du uj d+1 ~ 
N -(d+2)/{d+i) for fixed BC In 2D using kinetic theory 
and localization theory we expect localized modes for ui ~ 
lo^ = (In N)~ 1/2 , ballistic modes for u ~ wf = N~ 1/3 and 
diffusive modes for ui^ ~ ui ~ cj^. The diffusive contri- 
bution to total current will scale as Jq ~ (lniV) _1 ' 2 iV _1 . 
As argued above, the ballistic contribution depends on 
BCs with J B ~ N-^ 3 for fixed BC and J B ~ A^ 2 / 3 for 
free BC. Hence, adding all the contributions, we conclude 
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Table 1: The table gives the AT dependence of J and sum- 
marizes the main results of the paper. The error bar for the 
numerically obtained exponent values is of the order ±0.02. 
N.B: The system sizes used may be far from asymptotic. 



that asymptotically: J ~ (mA/)" 1 / 2 ^ 1 for fixed BC and 
J ~ N~ 2 / 3 for free BC. In 31? we expect that uj^ is in- 
dependent of N and states with to > lo^ are localized. 



jj? = N X I A are ballistic while 



Extended modes with < 
those with wf ~ ui ~ w^ are diffusive. The contribution 
from the diffusive modes scales as Jb ~ A^ 1 while the 
ballistic contribution again depends on boundary condi- 
tions with J B ~ iV~ 5 / 4 for fixed BC and J B ~ A^ 3 / 4 for 
free BC. We conclude that asymptotically: J ~ A^ 1 for 
fixed BC and J - A^ 3 / 4 for free BC. 

Pinned lattice: There is now a gap in the spectrum, 
starting from uj — 0, and thus there are no low frequency 
ballistic modes. In 2D localization theory then tells us 
that for sufficiently large N, all non-zero frequency modes 
become localized, hence we should get an insulator. In 3D 
there is a finite band of diffusive states and therefore we 
expect a normal conductor satisfying Fourier's law. 

Numerical results: We now check the above predic- 
tions through extensive numerical calculations. We mea- 
sure force-constants in units of k, masses in units of the 
average mass m, time in units of the inverse frequency 
fl^ 1 = (m/fc) 1 / 2 , displacements in units of the lattice 
spacing a, friction constant 7 in units of fhfl, and tem- 
perature in units of fha 2 £l 2 /ks- We set k = l,fn = 1 and 
we fixed AT = 1,7 = 1. For fixed BCs and pinned cases 
k' = 1. Different values of the mass variance A and the 
on-site spring constant k were studied for 2D and 3D lat- 
tices of different sizes 
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Here we report results for cases 
with the strongest disorder. The transmission coefficient 
TJv(w) can be expressed in terms of phonon Green's func- 
tions and can be accurately determined numerically using 
a transfer matrix representation [24] . By performing a 
discrete sum over the transmitting range of frequencies 
we evaluate the integration in Eq. (pi) to obtain J . For 
A^ > 64 in 3D the transfer matrix method has numerical 
problems and in those cases we performed nonequilibrium 
simulations to find the heat current. 

We also studied the properties of the normal modes of 
the isolated harmonic lattices with fixed BCs. We measure 
the degree of localization of a given mode by the inverse 

,2\2 



participation ratio defined as P 



En</(E r 



For an extended state P^ 1 is of order (D(N~ d ) while for a 
localized state, it is 0(1). 

In Fig. (fTl) we show the A^-dependence of the disorder 
averaged current J = \J\ in d = 2, 3. Our findings and the 
comparison with the heuristic predictions are summarized 
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Fig. 2: (i) Plot of T{ui). (ii) Plot of NT{u). The range of frequencies for which T(w) ~ 1/7V is indicated by arrowed lines and 
corresponds to diffusive modes . (iii) Plot of p{u)) for binary mass ordered lattice and single disordered sample, (iv) Plot of 
jyd-i x p-i £ Qr s j n gi e samples. Non-collapse of plots for different N in (iv) indicate localized states. 



in Table (flj . The results for the free BC case are consistent 
with the heuristic predictions while those for fixed BC 
show significant deviations. For the pinned system we 
obtain, in agreement with the theoretical prediction. In 
2D we find a strong decay of the current with N with 
/i = 3.7 suggesting that at larger system sizes we will get 
an exponential decay corresponding to a heat insulator. In 
3D we find \i = 1 implying that the system is a normal heat 
conductor. For the binary mass distribution we do not 
find a transition to insulating behaviour with increasing 
disorder. 

The iV-dependence of the disorder-averaged phonon 
transmission coefficient T(uj) = [T/v(w)]/A d_1 sheds addi- 
tional light on the nature of phonons in different frequency 
regions. For the case of fixed BCs we show in Fig. (pi) 
results for transmission, density of states and inverse par- 
ticipation ratios and from these we can see the range of 
allowed modes and their degree of localization. By plot- 
ting NT(lu) we identify the diffusive regime. From plots 
(iii-iv) we see that in both 2D and 3-D there are effectively 
two phonon bands, a remnant of the ordered binary mass 
case. In 2D the upper band is fully localized while in 31? 



there is a small number of localized states near the band 
edges. In 2D the lower band has extended states below a 
cut-off u>c which decreases slowly with N. In both 2D and 
3D the lower band has diffusive and ballistic states and 
the crossover scale uj^ decreases with system size. The ex- 
pected N— dependence of w^ and ui^ are difficult to verify 
at these system sizes. 

For the case of free BCs, we find that the values of T{uj) 
in the diffusive regime matches with those for fixed BCs 
but are completely different in the ballistic regime. This is 
seen in Fig. (pi) where we plot the effective mean free path 
ieff(w) = iVT(u;)/w d_1 in the low-frequency region. The 
difference between free and fixed BCs is larger in 2D than 
in 3D and explains the similar feature observed for the J 
values in Fig. (fTl). For free BC, Z c ff is roughly consistent 
with the kinetic theory prediction l~^ ~ N^ 1 + Ik~ (w) 
but the behaviour for fixed BC is very different. The inset 
of Fig. (pi) plots l e g for the equal mass ordered case and 
we find that in the ballistic regime it is very close to the 
disordered case, an input that we used in the heuristic 
derivation. The oscillations in the transmission for fixed 
BC arise from scattering and interference of waves at the 
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(b) 3D lattice 

Fig. 3: Plot of the effective mean-free path l c g — NT{uj)/oj d ~ 1 
in (a) 2D and (b) 3D with A = 0.8. The insets show £ cff for 
the ordered lattices with a single mass. An u~^ d+1 > behaviour 
is observed in a small part of the diffusive region. The fixed BC 
data is highly oscillatory and has been smoothed. 



interfaces. For the fully pinned case we find T(u>) 
in 2D and as T(u) ~ 1/N in 3D (24 . 



-bN 



In summary the analytic arguments show that the con- 
tribution of ballistic modes to conduction is dependent on 
BCs and is strongly suppressed for fixed BCs, the more 
realistic case. In 3D this leads to diffusive modes domi- 
nating for large system sizes and Fourier's law is satisfied. 
Thus a finite heat conductivity may be obtained without 
invoking anharmonicity as is usually believed to be neces- 
Our numerical results verify the predictions for 
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sary 

free BCs and we believe that much larger system sizes are 
necessary to verify the fixed BC results ( this is also the 
case in ID (18]). It is also possible that the assumptions 
made to predict the asymptotic N dependence of current 
is not correct. We are uncertain at present whether the 
disagreements are due to finite size effects or to the inade- 
quacy of the theory. A number of recent experiments have 
directly observed localization of electrons [I] , photons (2J , 
matter waves [5] and acoustic waves J3J. Our study shows 
clearly that the ballistic and diffusive contributions to the 
heat current would make it difficult to observe localization 
effects in heat conduction studies [6|[l4] . Our study of the 
pinned system gives us a clearer understanding of the role 
of low frequency modes in giving rise to diverging ther- 
mal conductivity and provides the simplest example of a 



3D deterministic system with diffusive transport. The 2D 
pinned system is a heat insulator and thin films attached 
on insulating substrates may show such behaviour. 
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